/*It ranks the combination of exclusion ranges and polynomial order by MSE*/
global data_in=20
global data_out=1180  
global bin=10
global binname=$bin 
global zoom_in_c300=80        
global zoom_out_c300=470
global range_c300   "incbin${binname}>=float(${zoom_in_c300}) & incbin${binname}<=float(${zoom_out_c300})" 

global low_range=-70
global high_range=20  
matrix MSE=J(4000,5,.)

global rr=1
/* Estimating Counterfactual - polynomial order*/			
foreach degree of numlist 3 (1) 7 {
	 global degree=`degree'

/* Estimating Counterfactual - exclusion region*/	
foreach low of numlist  ${low_range} ($bin) 0             {	
foreach high of numlist 0            ($bin) ${high_range} {
        global low=`low'
	    global high=`high'
	    matrix MSE[${rr},1]=${rr}
		matrix MSE[${rr},2]=${degree}
		matrix MSE[${rr},3]=${low}
		matrix MSE[${rr},4]=${high}
				
foreach num of numlist 1(1) 5 {
        
		use "${density_mse_dir}\density_${firm_name}_${name}_${year_name}_`num'_bin${binname}_in${data_in}_out${data_out}.dta", clear
		rename numbin${binname}  y_actual
	    qui  merge m:m  incbin${binname} using  "${density_mse_dir}\density_${firm_name}_${name}_${year_name}_exclude_`num'_bin${binname}_in${data_in}_out${data_out}.dta"			
		qui keep if $range_c300
		//use _exclude_1 to prdict _1
	    
		global notch_s=${cutoff}/${bin}  
		global low_s=${low}/${bin}    		
		global high_s=${high}/${bin}
		local y_actual="y_actual"	
		local incbin_name="incbin${binname}"
		local freq_name="numbin${binname}"

		gen round50=0
		replace round50=1 if mod(incbin${binname},50)==0

        gen Z100=0
        replace Z100=1    if mod(incbin${binname},100)==0   
			
	/* MATA ESTIMATION */		
		mata:mata clear
		mata:incbinname=st_local("incbin_name")   
		mata:freqname  =st_local("freq_name")  
		mata:yactual   =st_local("y_actual")
		mata:notch_s   =strtoreal(st_global("notch_s"))  
		mata:low=${low_s}
		mata:high=${high_s}
		mata:degree=${degree}
		mata:sf_s=$bin 
		
		
		* CALLING MATA PROGRAM TO ESTIMATE THE COUNTERFACTUAL *
	 if "${density}" =="DP_q"	& "${firm}"=="domestic" { 
		 mata:notch_mse_roundq(incbinname,freqname,yactual,notch_s,low,high,degree,sf_s)
		 }
	/* POPULATING mse */
	   local mse_`num'=${mse}
	} 
	/* data loop ends */
   local mse =`mse_1'+`mse_2'+`mse_3'+`mse_4'+`mse_5'
	   matrix MSE[${rr},5]=`mse'
	   global rr=${rr}+1	
	}	
	} 
	}   
	/* Degree loop ends */
   
        svmat MSE
		keep MSE* 
		rename MSE1  rr
		rename MSE2  degree
		rename MSE3  low
		rename MSE4  high 
		rename MSE5  mse

        keep if mse~=.
        sort mse  degree      	


